﻿using System;
using System.Collections.Generic;

namespace Allegro.Mathlib
{
    public static class Gamma
    {
        private const int ITMAX = 100;
        private const double EPS = 3.0e-7;
        private const double FPMIN = double.MinValue;

        private static double[] cof =
            {
                76.18009172947146,
                -86.50532032941677,
                24.01409824083091,
                -1.231739572450155,
                0.1208650973866179e-2,
                -0.5395239384953e-5
            };

        static public double P(double a, double x)
        {
            double gln;
            if (x < 0.0 || a <= 0.0)
            {
                throw new ArgumentException("Invalid arguments in function Gamma.P");
            }
            if (x < (a + 1.0))
            {
                double gamser;
                gser(a, x, out gamser, out gln);
                return gamser;
            }
            else
            {
                double gammcf;
                gcf(a, x, out gammcf, out gln);
                return 1.0 - gammcf;
            }
        }
        static public double Q(double a, double x)
        {
            double gln;
            if (x < 0.0 || a <= 0.0)
            {
                throw new ArgumentException("Invalid arguments in function Gamma.Q");
            }
            if (x < (a + 1.0))
            {
                double gamser;
                gser(a, x, out gamser, out gln);
                return 1.0 - gamser;
            }
            else
            {
                double gammcf;
                gcf(a, x, out gammcf, out gln);
                return gammcf;
            }
        }
        static void gser(double a, double x, out double gamser, out double gln)
        {
            double sum, del, ap;
            gln = gammln(a);
            if (x <= 0.0)
            {
                if (x < 0.0)
                {
                    throw new ApplicationException("x less than 0 in routine gser");
                }
                gamser = 0.0;
                return;
            }
            else
            {
                ap = a;
                del = sum = 1.0 / a;
                for (int n = 1; n <= ITMAX; n++)
                {
                    ++ap;
                    del *= x / ap;
                    sum += del;
                    if (Math.Abs(del) < Math.Abs(sum) * EPS)
                    {
                        gamser = sum * exp(-x + a * Math.Log(x) - gln);
                        return;
                    }
                }
            }
            throw new ApplicationException("a too large, ITMAX too small in routine gser");
        }
        static void gcf(double a, double x, out double gammcf, out double gln)
        {
            double an, b, c, d, del, h;

            gln = gammln(a);
            b = x + 1.0 - a;
            c = double.MaxValue;
            d = 1.0 / b;
            h = d;
            for (int i = 1; i <= ITMAX; i++)
            {
                an = -i * (i - a);
                b += 2.0;
                d = an * d + b;
                if (Math.Abs(d) < FPMIN)
                    d = FPMIN;
                c = b + an / c;
                if (Math.Abs(c) < FPMIN)
                    c = FPMIN;
                d = 1.0 / d;
                del = d * c;
                h *= del;
                if (Math.Abs(del - 1.0) < EPS)
                    break;
            }
            if (ITMAX > ITMAX)
            {
                throw new ApplicationException("a too large, ITMAX too small in gcf");
            }
            gammcf = exp(-x + a * Math.Log(x) - gln) * h;
        }

        static double gammln(double xx)
        {
            double x, y, tmp, ser;

            y = x = xx;
            tmp = x + 5.5;
            tmp -= (x + 0.5) * Math.Log(tmp);
            ser = 1.000000000190015;
            for (int j = 0; j <= 5; j++)
            {
                ser += cof[j] / ++y;
            }
            return -tmp + Math.Log(2.5066282746310005 * ser / x);
        }

        static double exp(double x)
        {
            double y = Math.Exp(x);
            return y;
        }
    }
}
